Modification of roton instability due to the presence of a second dipolar Bose-Einstein 
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We study the behavior of two coupled purely dipolar Bose-Einstein condensates, each located in 
a cylindrically symmetric pancake-shaped external confining potential, as the separation b between 
the traps along the tight confining direction is varied. The solutions of the coupled Gross-Pitaevskii 
and Bogoliubov-de Gennes equations, which account for the full dynamics, show that the system 
behavior is modified by the presence of the second dipolar BEC. For sufficiently small 6, the presence 
of the second dipolar BEC destabilizes the system dramatically. In this regime, the coupled system 
collapses through a mode that is notably different from the radial roton mode that induces the 
collapse of the uncoupled system. Finally, we comment on the shortcomings of an approach that 
neglects the dynamics in the z-direction, which is assumed to be a good approximation for highly 
pancake-shaped dipolar BECs in the literature. 

PACS numbers: 



I. INTRODUCTION 

Dipolc-dipole interactions are long-range and 
anisotropic and dominate the behavior of many liquids 
and solids such as ferrofluids and superfluid 3 He 
The condensation of 52 Cr atoms @, which have a large 
magnetic dipole moment compared to alkali atoms, 
paved the way for studying the physics of long-range 
interactions in a clean model system [l], While 
the dynamics of dipolar Bose-Einstein condensates 
(BECs) is, in general, governed by an interplay between 
the short-range s-wave interactions and the long-range 
dipole-dipole interactions, the s-wave scattering length 
can be tuned to vanish through the application of 
an external magnetic field in the vicinity of a Fano- 
Feshbach resonance @. This feature allows for the 
experimental realization of purely dipolar BECs. The 
present work investigates the behavior of two coupled 
dipolar BECs in a double-well type set up within the 
mean-field framework, which is expected to describe 
the key features of dipolar gases such as Cr BECs 
properly but not necessarily those of molecular samples 
such as RbK 0-0, 03 ■ Previous mean- field studies of 
single dipolar BECs in a pancake-shaped external trap 
predicted interesting features such as a red blood cell 
type shaped ground state density as well as collapse 
induced by radial and angular roton modes [ITrll3| . 

The behavior of dipolar BECs is even richer when a 
double well geometry is considered [T3 - tT6| . Describing 
the condensate by a single mean-field wave function, the 
existence of an instability island immersed in an other- 
wise stable region has been predicted to exist for certain 
parameter combinations [17]. Furthermore, macroscopic 
quantum self trapping, a phenomenon intensely studied 
for s-wave interacting BECs [l8l - l20j . has been predicted 
to occur for a dip olar BEC in a cigar shaped double- 
well potential jl4 [TEj . The transition from the macro- 
scopic quantum self-trapping to the Josephson oscilla- 
tion regime has been interpreted using a single two-mode 



model that treats the left well and the right well as be- 
ing occupied by macroscopic wave functions 'Fi and $2, 
respectively. Extensions to triple-well potentials, which 
provide a simplifying model of an optical lattice system, 
have also been considered [21[. Here, we model a two- 



well dipolar system, for which tunneling is assumed to be 
negligible, and solve a set of two coupled Gross-Pitaevskii 
(GP) and Bogoliubov-de Gennes (BdG) equations. Un- 
like the two-mode model eluded to above and unlike re- 
lated earlier studies [124131, our approach accounts for 
the full system dynamics within the mean-field frame- 
work. The dipoles in the two traps are assumed to be 
aligned along the tight confining direction and the system 
behavior is investigated as a function of the separation 
between the two clouds. 

The remainder of the paper is organized as follows. 
Sectionllll introduces the stationary and dynamical mean- 
field description of two coupled dipolar BECs. Section Hnl 
presents and interprets our numerical results of the sta- 
bility of the system as functions of the dipole strength, 
the aspect ratio and the separation between the two 
clouds. Lastly, Sec. II VI concludes. 



II. MEAN-FIELD DESCRIPTION 

A. Coupled Gross-Pitaevskii equations 

We consider two dipolar systems, each confined by a 
cylindrically symmetric external trap Vtj(p, z), j = 1 and 
2, 

V tj (p, z) = \mu 2 p (p 2 + A 2 (z - z jo f ) (1) 

with aspect ratio A, where A = u z /u pi and lo p and w z 
denote the angular trapping frequencies along the p and 
z directions, respectively. Here, we use cylindrical co- 
ordinates [r —{p,4>,z)\. In Eq. ([T]), m denotes the mass 
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of the dipoles and the Zjo denote the trap centers along 
the z-direction. We consider the same type of atomic 
species in both traps, e.g., 52 Cr in the same internal state. 
Throughout, we assume that the number Nj of dipoles 
in the j th trap (j = 1 or 2) is fixed, i.e., we assume that 
tunneling between the traps is absent. Our main interest 
is in determining the system behavior for aspect ratios of 
the order of 10 as the distance b, b — \z\q — Z20I, is var- 
ied. The parameter b determines the effective coupling 
between the two dipolar BECs. When b is infinitely large, 
the effective coupling vanishes and the system behaves 
like two independent dipolar BECs. When b is small, the 
effective coupling is strong and the system behavior is 
changed due to the long-range and anisotropic dipolc- 
dipole interaction between the dipoles located in the 
two traps. Throughout, we assume that the dipoles are 
aligned along the z-direction, so that the dipole-dipole in- 
teraction potential is given by Vdd(r) = rf 2 1 ~ 3 ™ s - , and 
s-wave interactions are neglected. Here, d is the dipole 
strength, r is the distance vector between the two dipoles, 
r = |r|, and •& is the angle between the z-axis and f. 

In the mean-field approximation, the two dipolar BECs 
are described by two coupled time-dependent GP equa- 
tions IH HI 



ih 



dt 



(Nj-1) / d 3 r'\* j (r' > t)\ 2 V dd (r-r') + 



N t / dV'|^(rV)| 2 y dd (r-r') 



(2) 



where (j, /) = (1, 2) and (2, 1), and H® denotes the single 
particle Hamiltonian, 
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V 2 + F tJ (p,z). 



(3) 



The coupling between the wave functions ^i(r, t) and 
1 t>2(r,t) arises due to the dipole-dipole interaction be- 
tween the two clouds, which is accounted for by the third 
term in the square bracket on the right hand side of 
Eq. The coupled mean field equations depend on 

four parameters, the aspect ratio A, the separation 6, and 
the dimensionless dipole strengths Di and D 2 , where 



Da 



d 2 {N 3 1) 



(4) 



Here, E p and a p denote the oscillator energy and length 
along the p-direction, E p = huj p and a p = ^Jh/(muj p ). 
The wave functions ^Sj are normalized according to 
J d 3 f\ 1 $j(r, t)\ 2 = 1. Since the confining potential and 
the dipole-dipole interaction potential are cylindrically 
symmetric, we can write the wave functions as ^ j{r, t) = 
^j(p,z,t)exp(ik(f>), where k = 0, ±1, ±2, • • •. 

The ground state solution of Eq. ^ can be obtained 
by solving the coupled time-independent GP equations 



self-consistently [26|. To this end, we set k = and 
write ^j(r,t) = ipjipi z ) exp(—ipjt/h), where the pj de- 
note the chemical potentials corresponding to the ground 
state solutions ip®(p,z). We solve the coupled time- 
independent GP equations self-consistently by evolving 
an initial state in imaginary time until convergence is 
reached [27| . Our numerical implementation exploits 
the cylindrical symmetry of the system and uses a two- 
dimensional grid in the p- and z-directions [IH, [28j . 

In addition to the stationary ground state wave func- 
tions ip® and the chemical potentials pj, we determine 
the total energy per particle E to t/N, 



a — i n J 



i=i,2 



£ V2 + ^' Z) + 



I d 3 f'\^(f')\ 2 V dd (f-f') 



d'r / d d r'|^(r")|V dd (f-r')|V 2 U (r)| 2 ,(5) 
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where N = Ni + Here, we used ip®(p,z) = "4>j(r) 
for notational convenience. The first, second, and third 
terms in the square bracket on the right side of Eq. ([5]) 
give rise to the kinetic energy per particle Eki n /N, the 
trap energy per particle E trap /N ', and the on-site dipole- 
dipole interaction energy per particle Edd,on/N , respec- 
tively. The off-site dipole-dipole interaction energy per 
particle Edd,off/N is given by the last term on the right 
hand side of Eq. ([5]) . 

We note that the system considered here can be 
viewed as a variant of the first mean-field study of two- 
component dipolar BECs by Goral et al. [26| who consid- 
ered the limiting case of vanishing separation and spher- 
ical confinement. However, as opposed to two dipolar 
BECs aligned along the same direction as considered 
here, they considered two oppositely polarized BECs. 
We have checked for selected cases that our solutions for 
the coupled stationary GP equations agree with those 
reported by Goral et al. 



B. Coupled Bogoliubov de Gennes equations 



To analyze the dynamical behavior of the system, we 
write [H [33 

*j (r, t) = [^{r) + 5il>j (r, t)] exp(-»^i/ft), (6) 

where j = 1 and 2, and the S^j(r,t) denote the pertur- 
bation of the dipolar BEC located in trap j. Following 
the literature [301 ] , we write the perturbations in terms of 
the Bogoliubov particle and hole excitations Uj and Vj , 

5vpj(r, t) — Uj(r) exp(— iujt) + v*(r) exp(zwi). (7) 

Plugging Eqs. ^ and (J7J into Eq. © and keeping terms 
up to the first order in Uj and Vj, we find, after equating 
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the coefficients of exp(iwi) and exp(— iujt), a set of two 
coupled BdG equations, 



h 2 LU 2 f 3 (r)=A J {r)A J {r)f J {r) + 
2{N j - 1) I #?'ilP j {i f ')f i (?')V4 A <? ~ f) + 

d 3 f'^{f')fi(r')V dd (r-f) 



2Ni 



$(0. (8) 



where, as before, (J, I) = (1,2) and (2,1). In deriv- 
ing Eq. ©, we assumed that the tpj(r) are real. The 
functions fj{r), fj(r) = Uj{r) +Vj(r), represent the den- 
sity perturbation for the dipolar BEC located in the j th 
trap [28|, [3l| . This becomes clear if we calculate the den- 
sity using Eq. ©. Assuming that the Uj and Vj are 
real and keeping only the lowest order correction, we ob- 
tain | tyj (r, t) | m \ipj(r)\ +2cos(u>t)ipj(r)fj(r). Because 
of the cylindrical symmetry of the problem, the density 
perturbations or eigen modes fj(r) can be written as 
fj(p,z)exp(ik(f>), k — 0,±1,±2, • • • [13]. The operators 
Aj (r) in Eq. © operate on everything to their right and 
are given by 



Mr) = H°(r) + (N S -I) J dV'|^(r")| 2 V dd (f - r ') + 

N t [ d 3 r"|^ (f')|V d d(f-f')-^-(9) 



We solve Eq. ([8]) for the eigen frequencies uj and the 
density perturbations fa and fa for various k, k = — 4, 
using the Arnoldi method [33f. Our implementation fol- 
lows that discussed in Ref. [28| for a single component 
dipolar BEC. In particular, we construct a vector from 
the density perturbations fa and fa, and then proceed as 
in the single component case. The excitation frequencies 
uj allow for the determination of the dynamical stability 
of the system. A positive uj 2 , or real uj, signals that the 
system is dynamically stable with respect to the associ- 
ated density oscillation. A negative uj 2 , or imaginary uj, 
in contrast, signals that the system is dynamically unsta- 
ble with respect to the associated density oscillation. As 
detailed further in Sec. IIII1 the dynamical instability of 
the two well dipolar system can, depending on the sys- 
tem parameters, be triggered by either a k = mode or 
a finite k mode. 
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FIG. 1: (Color online) D-versus-l/& phase diagram for (a) 
A = 7, (b) A = 8, (c) A = 9, and (d) A = 15. The solid lines 
separate the dynamically stable region from the dynamically 
unstable region (labeled as "U"). Circles, squares and dia- 
monds indicate that the system collapses through the soften- 
ing of a k = 0, k — 2 and k — 3 excitation mode, respectively. 
In panels (a), (b), and (c), the dashed lines separate the me- 
chanically stable region from the mechanically unstable region 
while the dash-dotted lines mark the boundary between re- 
gions where density maxima are located at p = (labeled 
as "So") and where density maxima are located away from 
p = (labeled as "S>o"). For A = 8, panel (b), dash-dash- 
dotted lines indicate, from bottom to top, 10%, 20%, 30%, 
and 40% difference between n,(0) and rij(pmax) (see text). 
The oscillator length a z is defined through a z — ^/fi/ (mu> z ). 



III. RESULTS 

As discussed in Sec. [Til the coupled GP equations de- 
pend on four parameters. To reduce the parameter space, 
we set Ni = N2 = N/2, and investigate the system prop- 
erties as a function of A, b, and D (D = D\ = Z^)- 
Figures H{a)-(d) show the D- versus- 1/6 phase diagram 
for A = 7, 8, 9 and 15, respectively. The solid lines sepa- 
rate the dynamically stable region from the dynamically 
unstable region. The symbols indicate the mode through 
which the system becomes unstable; circles, squares and 



diamonds correspond to the k = 0, 2 and 3 mode, re- 
spectively. For a fixed separation b, the dynamically 
stable region increases with increasing A. This behav- 
ior is well known for a single dipolar BEC in a pancake- 
shaped trap [25;, 34]. As the aspect ratio A increases, the 
dipole-dipole interaction becomes effectively more repul- 
sive. Figure [JJ shows that the dynamically stable region 
decreases with decreasing separation, i.e., increasing 1/6, 
for fixed A. This decrease of stability is attributed to 
the presence of the second cloud. Since the dipoles are 
aligned along the z-axis, the dipole-dipole interaction be- 
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FIG. 2: (Color online) Asymmetry parameter B as a function 
of 1/6 for A = 7 (solid line), A = 8 (dashed line), A = 9 (dot- 
dot-dashed line), and A = 15 (dash-dotted line). For a given 
A and b, the D value is chosen to be approximately equal to 
0.9D cr . 



tween the two neighboring clouds, Edd,off, becomes more 
attractive as the separation between the clouds is de- 
creased. 

The critical dipole strength D cr , defined as the D 
value for which the system becomes dynamically un- 
stable, changes in an interesting manner with increas- 
ing a z /b. For small a z /b, D cr varies slowly. Around 
a z/b > 0.11,0.13,0.14, and 0.17, D cr varies compara- 
tively fast for A = 7, 8, 9, and 15, respectively. Finally, 
for larger a z /b, D cr varies again comparatively slowly. In 
order to understand this dependence of D cr on a z /b, we 
analyze the ground state density of the system. 

Dashed lines in Figs. [Tfa)-(c) show the mechanical in- 
stability line. The mechanical and dynamical instabil- 
ity lines nearly coincide when the ground state densi- 
ties 1^1 2 are approximately Gaussian shaped (labeled as 
"So") but deviate when the densities iV'jl 2 have a so- 
called red blood cell type shape (labeled as "S>o"), i.e., 
when the density maxima are located at p > 0. The dash- 
dotted lines separate the two types of densities, which are 
determined by analyzing the integrated densities rij(p), 
where n.j(p) — 2ir J (r)\ 2 dz and j = 1 and 2, along 
the p direction. If n,j(0)/nj(p max ) < 0.98, where p max is 
the p value at which the integrated density rij (p) has its 
maximum, then we call the density red blood cell shape; 
otherwise we call it Gaussian. Figures Q^a) and (b) show 
that the ground state density near the instability line has 
red blood cell type structure for fairly large separation. 
When the inverse separation a z /b has increased to about 
0.1 — 0.13, the red blood cell type structure disappears. 
D cr changes more rapidly for intermediate a z /b values 
when the ground state density is Gaussian. Figures [Tf a) 
and (b) suggest that the deformation of the ground state 
density away from the simple Gaussian like profile leads 
to a significant stabilization of the system. This interpre- 
tation is supported by the fact that the deformation, or 
the red blood cell type structure, becomes comparatively 
more pronounced as a z /b increases from to about 0.13, 
as indicated by the dash-dash-dotted lines in Fig. [TJb) . 




FIG. 3: (Color online) (a) Energy contributions Ei/N as 
a function of D for A = 7 and a z /b = 1 / (5%/A) « 0.076 
(010 w — 6.61a z and Z20 ~ 6.61a z ). The solid, dashed, 
dash-dotted, dotted, and dash-dash-dotted lines show Etot/N, 
Ekin/N, Etrap/N, E dd , on j TV , and E dd , off /N, respectively, 
(b) The energy contribution E d d,off/N is shown as a function 
of D for A = 7 and a z /b = (solid line), a z /b = 1/{WVX) ^ 
0.038 (dashed line), a z /b = 1/(5 V / A) ~ 0.076 (dash-dash- 
dotted line), a z /b = l/(3.33\/A) « 0.113 (dotted line), and 
a z /b — 2/(5\/A) ~ 0.151 (dash-dotted line), respectively. 



For larger aspect ratios [see Fig. [IJc)], the ground state 
density in the dynamically stable region has red blood cell 
type structure only in a tiny region in the vicinity of the 
instability line around a z /b « 0.13. Even so, for slightly 
larger a z /b values, D cr varies more rapidly. For A = 15 
[see Fig. H{d)], the red blood cell type structure exists 
in an even smaller region around {D,a z /b) w (90,0.17), 
where the system collapses through a k — 2 mode. Al- 
though much less pronounced, D cr varies more rapidly 
for slightly larger a z /b values. Interestingly, the density 
deviates again from the simple Gaussian type shape in a 
tiny region around (D,a z /b) w (35,0.235) 

To further characterize the ground state density, we 
define the quantity Bj, which measures the asymmetry 
of the density of the cloud located in the j th trap, i.e., 
the density asymmetry about the trap center, 



Bj = 2tt 



pdpdz\^(r)\< 



OO POO 



pdpdz\^(f}\' 



(10) 



where j = 1 and 2. For equal number of dipoles in both 
traps, as considered in this paper, one has B\ = —B-2 
and we define B = \Bj\. When the separation b is 
large, the ground state densities \ip® | 2 are symmetric and 
B = 0. However, as b decreases, the interaction be- 
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tween the dipoles located in the two traps leads to an 
increased density between the trap centers and thus to 
finite B values. To quantify how B changes with decreas- 
ing b, we move along a "trajectory" in the D- versus- 1/6 
phase diagram for fixed A (see Fig. [IJ that lies 10% be- 
low the solid line, i.e., we choose D s» Q.9D cr for fixed 
A and b. Figure [5] shows that the asymmetry parameter 
B is, to within our numerical accuracy, identically zero 
for a z /b < 0.05 — 0.075 for the trajectories investigated. 
Figure [5] shows that it takes a certain critical attraction 
before the system breaks the symmetry of the ground 
state density. For large separations, the energy is mini- 
mized for densities symmetric about ZjQ. For smaller b, 
however, the off-site interaction is attractive enough to 
deform the ground state densities \ipj \ 2 . 

To gain further insight, Fig. (3{a) shows the energy per 
particle, E to t/N, as well as the individual energy contri- 
butions as a function of D for A = 7 and a z /b « 0.076. 
Etot/N (solid line), E tra p/N (dash-dotted line) and 
Edd,on/N (dotted line) increase with increasing D while 
Ekin/N (dashed line) and Edd,off/N (dash-dash-dotted 
line) decrease with increasing D. Since the dipoles are 
aligned along the z-direction and the system is pancake- 
shaped, Edd,on is effectively repulsive. To minimize 
Edd,on, the dipoles try to spread out, which reduces Eki n 
but increases Et ra p as D increases. Edd,off decreases 
(i.e., becomes more negative) with increasing D since the 
attraction between the dipoles located in the two traps 
increases. This qualitative behavior remains the same as 
the separation b between the traps decreases. For fixed 
D and A, the energy contribution that changes the most 
is Edd,off- As shown in Fig. [3jb), Edd,of}/N decreases 
appreciably as the separation b decreases, which can be 
attributed to the attraction between dipoles located in 
the two different clouds. The fact that Edd.off/N is the 
energy contribution that changes the most as b decreases 
emphasizes that the decreased stability is driven by the 
dipole-dipole interactions between the two clouds. 

We now investigate how the collapse mechanism 
changes with b. For infinite separation, each BdG excita- 
tion frequency is doubly degenerate and the system col- 
lapses through a k — mode when the ground state den- 
sity has Gaussian shape and through a k > mode when 
the ground state density has red blood cell shape fill ]. 
Figured] shows the seven lowest BdG eigen frequencies as 
a function of D for k = 0, A = 15, and a z /b = 0. In this 
case, the system collapse is triggered by a radial roton 
mode. Figure [5] shows the eigen mode f± corresponding 
to the lowest BdG eigen frequency for D « 144 (see tri- 
angle in Fig.@|. The density oscillation /i has six nodal 
lines that are separated by approximately 0.8a p . This 
nodal line spacing agrees quite well with A p /2, where 
A p is the wavelength expected for a radial roton mode, 
A p « 27ra z « 1.62a p [Hill. 

As b decreases, the collapse mechanism of the system 
changes due to the attractive dipole-dipole interaction 
between the two clouds. Figure |6] shows the fourteen low- 
est BdG eigen frequencies as a function of D for k — 0, 




FIG. 4: Solid lines show the seven lowest non-zero BdG exci- 
tation frequencies uj as a function of D for A = 15, a z /b = 0, 
and k — 0. The triangle marks the BdG eigen frequency 
for which Figs. [5] and [8] show the corresponding eigen mode 
fi (p, z) and ground state density, respectively. 




FIG. 5: (Color online) Contour plot of the BdG eigen mode 
fi(p,z) for A = 15, a z /b = 0, D « 144, and fe = 0. The 
contours are equally spaced. Solid and dashed lines indicate 
positive and negative /i (p, z) values, respectively. The dash- 
dotted lines indicate the nodal lines and the dotted line marks 
the trap center. 




% 10 20 30 40 50 
D 

FIG. 6: (Color online) Solid and dashed lines show the 
in-phase and out-of-phase BdG excitation frequencies tj as 
a function of D for A = 15, a z /b = 4/(5^) « 0.21 
(210 ~ — 2.42a z and 220 ~ 2.42o z ), and k — 0. The triangle 
and square mark the BdG eigen frequencies for which Fig. [7] 
shows the corresponding eigen modes /i(p, z) and /2(p, z)\ the 
corresponding ground state density is shown in Fig. [8] 
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FIG. 7: (Color online) Contour plot of (a) and (b) the bi- 
phase BdG eigen modes f\(p,z) and fi(p,z), respectively, 
and (c) and (d) the out-of-phase BdG modes /i (p, z) and 
f'2(p,z), respectively, for A = 15, a z /b = 4/(5v / A) ~ 0.21 
(zio « -2.42a z and z 20 « 2.42a*), D w 52.02, and fe = 0. 
The contours are equally spaced. Solid and dashed lines in- 
dicate positive and negative fi(p,z) and f2(p,z) values, re- 
spectively. The dash-dotted lines indicate the nodal lines and 
the dotted lines mark the trap center. The eigen frequency 
corresponding to panels (a) and (b) is marked by a triangle 
in Fig. [6] while that corresponding to panels (c) and (d) is 
marked by a square. 



A = 15, and a z /b « 0.21. Solid and dashed lines in- 
dicate that the BdG eigen frequencies correspond to in- 
phase and out-of-phase density oscillations of the BECs 
located in the two traps. For D = 0, the eigen spectrum 
consists of degenerate pairs. As D increases, the in-phase 
and out-of-phase frequency pairs decouple. This behav- 
ior is analogous to that of a symmetric one-dimensional 
double-well potential. In the weak coupling regime (high 
barrier), the tunneling splitting is small and the eigen 



spectrum consists of nearly degenerate pairs. As the 
coupling increases, the eigen frequency pairs decouple. 
Figure [6] shows that the lowest in-phase eigen frequency 
approaches zero for D ss 52.9. FiguresJTJa) and (b) show 
the corresponding eigen modes for a slightly smaller D 
value (see triangle in Fig. |6] the corresponding ground 
state density is shown in Fig. |8]). As a result of the 
attractive off-site dipole-dipole interaction, the in-phase 
eigen modes /i [Fig. EJa)] and f% [Fig. Efb)] just prior to 
collapse are slightly asymmetric around the trap centers 
z\o and Z20, respectively. The in-phase eigen modes /1 
and /2 have three nodal lines, whose separation increases 
slightly with increasing p. This suggests that the coupled 
dipolar BEC system does not, like the uncoupled system 
(see Figs. |4]and [5}, collapse through a "pure" radial ro- 
ton mode. The radial roton mode can be interpreted as 
being the result of the formation of a pattern along the 
p-direction whose size is governed by a z . In the presence 
of the second dipolar BEC, the separation b between the 
two clouds sets another length scale. For the parameters 
in Figs. |6] and we have b = 1.25a p « 4.84a z . Thus, the 
dynamics along the p-direction is governed by an inter- 
play of the length scales a z and 6, resulting in modes /1 
and ji that have neither the characteristic features of a 
"pure" radial roton mode nor those of a "pure" breathing 
mode of the entire two-cloud system. 

For comparison, Figs. [7tc) and (d) show the eigen 
modes f± and fa corresponding to the second lowest out- 
of-phase frequency for D « 52.02 (see square in Fig. [6]). 
For vanishing D, the corresponding eigen frequency is de- 
generate with the in-phase eigen frequency that, for finite 
D, triggers the collapse. The nodal pattern of the out-of- 
phase eigen modes [Figs.[7jc) and (d)] is distinctly differ- 
ent from that of the in-phase eigen modes [Figs.JTJa) and 
(b)], underlining the fact that the eigen frequency pairs 
decouple as the coupling between the clouds increases. 

Finally, we investigate the system behavior assuming 
that the dynamics in the z-direction is frozen, i.e., we 
write [H-|2| 

* j (r,t)=$ j (p,4>,t)4>° j (z), (11) 

where the (j>j(z) denote the one-dimensional harmonic 
oscillator ground state wave functions of the j th trap in 
the z-direction. For the parameter combinations inves- 
tigated, the ground state densities obtained using the 
frozen z-dynamics approach show, just as the densities 
obtained using the full mean-field dynamics, Gaussian 
and red blood cell type structures in the vicinity of the 
mechanical and dynamical instabilities. However, these 
structures appear at different (D, b) combinations for the 
two different approaches. We find that the variational 
wave function given in Eq. (|lll) predicts the mechani- 
cal instability to set in at much larger D values than 
predicted by the mean-field wave function that accounts 
for the full dynamics. Moreover, the frozen z-dynamics 
approach predicts a fairly smooth decrease of D cr with 
increasing 1/6 (for a z /b < 0.25) and does not reproduce 
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the relatively steep drop ofD cr around a z /b = 0.1 — 0.175 
discussed in the context of Fig.[T] Figures ISlTTUlexemplar- 
ily illustrate these findings. 

To start with, we analyze the energetics for A = 15, 
a z /b ss 0.21 and D » 52 (see triangle in Fig. [6]). We find 
that the ground state energy obtained for the variational 
wave function, Eq. (fTTj) . is about 5% higher than the 
exact mean-field energy. While the total energy agrees 
fairly well, the kinetic energy Eki n differs by about 40%, 
suggesting that the description based on the frozen z- 
dynamics is not flexible enough to describe all features 
of the system qualitatively correctly. Indeed, the frozen 
z-dynamics approach predicts the dynamical instability 
to occur at D cr fts 515 (see Fig. [HI, i.e., D cr predicted by 
the frozen z-dynamics is about ten times larger than D cr 
predicted by the full mean-field dynamics. 

Solid lines in Figs. HJa) and (b) compare the contour 
plots of the ground state density obtained accounting for 
the full dynamics and assuming frozen z-dynamics, re- 
spectively, for A = 15, a z /b 0.21 and D values that are 
slightly smaller than the respective D cr , i.e., D « 52 in 
Fig. Ha) and D « 512 in Fig. [S(b). The ground state 
density obtained assuming frozen z-dynamics [Fig. [SJb)] 
is significantly more extended in the p-direction and less 
extended in the z-direction than that obtained account- 
ing for the full dynamics [Fig. UK a)]. For comparison, 
dashed lines in Figs. Eta) and (b) show the ground state 
density for the same A value (i.e., A = 15), but 1/6 = 0. 
The D values are, as for the densities shown by solid lines, 
chosen to be slightly smaller than the respective D cr val- 
ues, i.e., D » 144 in Fig.^a) and D « 608 in Fig.^b). 
Compared to the densities for finite separation, those for 
infinite separation are more extended in the p direction. 

Figure [5] shows the BdG eigen spectrum obtained as- 
suming frozen dynamics in the z-direction as a function 
of D for k = 0, A = 15, and a z /b « 0.21. A compari- 
son of Figs. [9] and [6] shows that the spectrum obtained 
based on the frozen z-dynamics reproduces that obtained 
based on the full dynamics qualitatively but not quanti- 
tatively. In particular, the frozen z-dynamics approach 
predicts D cr w 515, compared to D cr w 52.9 obtained 
using the full dynamics approach. Figure [lU] shows the 
eigen mode /i for A = 15, a z /b w 0.21, D « 512 (which 
is just a bit smaller than D cr ), and k — for the lowest 
eigen frequency. The eigenmode possesses 13 nodal lines, 
which are approximately equally spaced (the spacing is 
about 0.55a p for the first 8 or 9 nodal lines and slightly 
larger for the last 5 or 4 nodal lines) , indicating that the 
collapse is triggered, according to the frozen z-dynamics 
approach, by a radial roton mode and not, as predicted 
by the full mean-field dynamics, by a mode that is nei- 
ther a pure radial roton mode nor a pure breathing mode 
[see discussion around Fig. Ufa)]. Although it might be 
expected intuitively that an aspect ratio of A = 15 is 
sufficiently large to treat the system as effectively one- 
dimensional, our analysis shows that this is not the case. 
Our study suggests that caution needs to be exercised 
when the dynamics of coupled pancake shaped traps is 
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FIG. 8: (Color online) Contour plot of the ground state den- 
sity for A = 15 calculated (a) accounting for the full dynamics 
in the z-direction and (b) assuming frozen dynamics in the z- 
direction. The solid lines correspond to a z /b m 0.21 and (a) 
D ~ 52 and (b) D « 512. The dashed lines correspond to 
a z /b = and (a) D « 144 and (b) D w 608. The contours 
correspond to 90% (centermost contour), 70%, 50%, 30% and 
10% (outermost contour) of the peak density. The dotted 
lines mark the trap center. 



treated within a variational approach. Future studies 
need to extend the analysis to even higher A to make 
direct contact with Refs. EB-EmL 



IV. SUMMARY 

We studied the behavior of two coupled dipolar BECs, 
each located in a cylindrically symmetric external con- 
fining potential, as the separation b between the traps 
along the tight confining direction is varied. The dipoles 
are aligned along the z-direction and s-wave interactions 
are neglected. The number of dipoles in each trap is 
conserved separately, i.e., tunneling between the traps is 
neglected. The solutions of the coupled GP equations 
show that the system behavior is modified by the pres- 
ence of the second dipolar BEC. As the separation is 
decreased from infinitely large values to a value of about 
5 or Aa z , initially the collapse behavior changes little and 
then significantly below a certain critical separation. For 
separations smaller than this critical separation, the pres- 
ence of the second dipolar cloud destabilizes the system 
dramatically compared to the case where the traps are 
infinitely far apart. For certain parameter combinations, 
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FIG. 9: (Color online) Solid and dashed lines show the in- 
phase and out-of-phase BdG excitation frequencies uj obtained 
assuming that the z-dynamics is frozen as a function of D 
for A = 15, a z /b m 0.21, and k — 0. The triangle marks 
the BdG eigen frequency for which Figs. [5] and QJj] show the 
corresponding density and eigen mode fi(p,z), respectively. 
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FIG. 10: (Color online) Contour plot of the BdG eigen mode 
fi{p,z) for A = 15, a z /b « 0.21, D « 512, and k = obtained 
assuming that the 2-dynamics is frozen. The contours are 
equally spaced. Solid and dashed lines indicate positive and 
negative fi (p, z) values, respectively. The dash-dotted lines 
indicate the nodal lines and the dotted line marks the trap 
center. The eigen frequency corresponding to this eigen mode 
is marked by a triangle in Fig. O 



we find that the so called red blood cell type density be- 
comes more pronounced or appears due to the presence 
of the second dipolar BEC. For infinitely large separa- 
tion, each BdG frequency is doubly degenerate. As the 
separation is decreased, the solutions of the coupled BdG 
equations show that the eigen frequency pairs decouple 
into two eigen frequencies corresponding to in-phase and 
out-of-phase density oscillations of the BECs located in 
the two traps. When the separation between the traps is 
large, the system collapses through a radial roton mode if 
the ground state density is Gaussian shape and through 
an angular roton mode if the ground state density is red 
blood cell shape, similar to the case of a single dipolar 
BEC [ll|. For relatively small separation, in contrast, the 
system collapses through a mode that is notably different 
from the radial roton mode that induces the collapse of a 
single dipolar BEC. For comparison, we also considered 
a simplified description in which the dynamics in the z- 
direction is assumed to be frozen. Compared to the full 
mean-field description, the simplified description, which 
is used frequently in the literature [23 - 124) . reproduces 
some features qualitatively but not quantitatively. We 
conclude that the frozen z-dynamics approach is inade- 
quate to quantitatively describe certain aspects of purely 
dipolar BECs, including the collapse, even if the aspect 
ratio is fairly large. 



During the final stage of preparing this manuscript for 
submission, we became aware of a related study by Wil- 
son and Bohn [3(| that considers the dynamics of an ar- 
ray of dipolar pancake-shaped BECs at various levels of 
approximation. 



Support by the NSF through grant PHY-0855332 is 
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